Parameter distributions of Keplerian orbits 
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O ' ABSTRACT 

a: 

i Starting with just the assumption of uniformly distributed orbital orientations, we 

derive expressions for the distributions of the Keplerian orbital elements as functions 
of arbitrary distributions of eccentricity and semi-major axis. We present methods for 
finding the probability density functions of the true anomaly, eccentric anomaly, orbital 
[ radius, and other parameters used in describing direct planetary observations. We 

also demonstrate the independence of the distribution of phase angle, which is highly 
significant in the study of direct searches, and present examples validating the derived 

q ' expressions. 

Subject headings: celestial mechanics — methods: analytical — methods: statistical — 
^ , planets and satellites: detection 

1. Introduction 



Brownl (J200J) introduced the concept of completeness to study the selection effects introduced 
by observatory architectures on direct searches for sub-stellar companions. Assuming distributions 
for semi-major axis and eccentricity of planetary orbits, Brown calculated the probability that 
a co mpanion would fall outside the telescope's central obscuration during an observation of a 



star. iBrownl (|2005l ) subsequently expanded this concept to include the selection effects due to 
the photometric restrictions on observability introduced by telescope optics, and IBrownl ((2009) 
demonstrated how completeness could be evaluated for indirect companion detection methods such 
as astrometry. Completeness has also been extended to account for multiple observations of one 
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star at different times (jBrown and Soummerll2010l ). and has been utilized in mission analysis an d 
development for a variety of proposed exoplanet observatories (jSavranskv et al.ll2010l ; lBrownl l2009). 



The direct detection (imaging) completeness is evaluated by assuming that a companion will 
be observable if its angular separation from the star is greater than the observatory's inner working 
angle (IWA), and illuminated such that the difference in brightness between star and companion 
(Amag) is below a threshold value, called the limiting Amag, or Amago. The IWA represents the 
minimum angular separation between the telescope center-line and detectable objects on the sky. 
It is determined by the size of a central obscuration, or the capability of adaptive optics systems to 
remove light from certain areas of the image plane, or the size and geometry of external occulting 
optics. Amago represents the point where systematic errors produce unresolvable confusion between 
planet signal and background noise. To calculate the completeness, probability distributions (or 
constant values) are assumed for planetary orbital elements and physical properties (see ^2]). A 
large, equal number of samples is generated from each distribution, and the star-planet angular 
separation and Amag are calculated for each set of samples. When binned in a two-dimensional 
histogram, these generate a density function representing the probability that a planet drawn at 
random from the assumed population will have a given angular separation and Amag. Integrating 
over this density yields a cumulative distribution function (CDF), which can be used to determine 
the probability that an observatory with a given Amago an d IWA, observing a specific star once, 
will be able to detect a planet belonging to the assumed population. 

The procedure described above is a Monte-Carlo sampling of a bivariate distribution function 
of non-independent arguments (since star-planet separation and Amag are functions of the same 
parameters). This means that to find any one point (or section) of the completeness distribution, 
it is necessary to sample it completely. Because of the relatively high dimensionality of the initial 
parameter space and wide range of values certain parameters can tak e , com plete sampling requires 
a large number of Monte-Carlo trials. The first simulation in iBrownl (120051 ). for example, includes 
100 million samples, and it can be shown that certain low probability areas of the function are 
under-sampled. 

Any alternate method of sampling completeness requires at least some knowledge of it s density 
function. In particular, Markov chain methods such as Metropolis-Hastings (|Hastingslll970l ) perform 
significantly better if the proposal distribution (a function used to 'propose' new samples that are 
then either accepted or rejected) closely approximates the target distribution. A special case of 
Metropolis-Hastings, known as Gibbs sampling, can be used to generate a sequence of samples from 
the joint distribution of two variables if their conditional probabilities are known. Our goal is to 
derive the distribution functions of the arguments to the completeness functions. 

We do so by starting with an ensemble of Keplerian orbits whose orientation is uniformly 
distributed in space, and deriving the distribution functions of these orbits' Keplerian parameters. 
Rather than constraining the population of orbits, we make no assumptions as to the distribution of 
orbital semi-major axis and eccentricity. This allows us to derive the completely general distribution 
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functions presented in £j3j Building upon this, we further consider parameters related to direct 
planetary observations in and make the discovery that the planetary phase angle (star-planet- 
observer angle) is independent of any of the orbital parameters. This makes it possible to write 
distribution functions for quantities directly related to the two parameters of the direct detection 
completeness function. It is important to note that, while the specific application explored here 
is direct imaging, the distributions derived in this paper are more broadly applicable to exoplanet 
studies in general. For example, Keplerian fits are often employed in doppler spectroscopy surveys, 
making these derivations useful for inferring the true distributions of orbital parameters derived 
from radial velocity data sets. Similarly, statistical analyses play an importa nt role in other met hods 
of exoplanet study, including transit photometry and microlensing surveys (jGould et al.ll2010l ). 



2. Modeling Keplerian Orbits 

In this section, we define the basic orbital parameters we will use throughout to describe 
companion orbits, and derive the relationships between these parameters and quantities that may 
be observed via direct planet searches. We begin our derivation with three basic assumptions: First, 
that planetary orbits may be approximated as closed Keplerian orbits with negative specific orbital 
energy. Second, that the distribution of orbital orientations with respect to an observer is uniform 
over a spherical volume. Finally, that any given star will have either one or no companions. The 
first two assumptions allow us to describe any orbit with the parameter set (a, e, ip, 6, 4>) where a is 
the semi-major axis, e is the eccentricity, and ip,0,4> are Euler angles determining the orientation 
of the orbit in the observer's reference frame. The position of the planet on its orbit at the time of 
observation is given by the true anomaly v. The third assumption allows us to describe orbits in 
an exact fashion (via the Kepler solutions to the two-body problem), but means that the derived 
equations will fail to capture effects due to planet-planet interactions such as resonant orbits. 

The observer's reference frame, X = (0,x,y,z), is defined via a Cartesian coordinate set with 
the origin O placed at the location of the observed star, and with the observer a distance d from 
O along the — z axis. A face-on orbit thus lies in the (x, y) plane, while an edge-on orbit lies in 
any plane orthogonal to the face-on one. Starting with an orbit lying in the (y, z) plane, with the 
planet's path counter-clockwise about the x axis, we apply a 3-1-3 rotation using the angle set 
(tp, 6, 4>). By the second assumption above, the angles ip and 4> are uniformly distributed in [0, 2n], 
while 9 is sinusoidally distributed in [0, 7r) — i.e., 6 is drawn from cos -1 ([/([0, 1])) where U is a 
uniform distribution. 
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The position of the planet with respect to its parent star, at the time of observation, is thus: 





COS 4> 


sin0 


o~ 




"l 










cos ip 


sin^ 


o" 







T p/s — 


— sm<p 


COS (j) 










cos 9 


sin 9 




— sini/' 


cosip 







rsini^ 










1 







— sin 9 


cos 9 










1 




r cos v 



(1) 



r (cos v sin 9 sin + sin v[cos 9 cos ?/> sin <j) + cos sin ip)) 
r(cos i/ cos 4> sin + sin z^(cos 9 cos cos ip — sin sin ip)) 
r(cos cos f — cos ip sin sin v) 

where r is the distance between the planet and the star, 

a(l-e 2 ) 
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Note that w e have used left-handed rotations here to conform to conventions in the literature (see, 
for example, Vinti (jl998l )). However, because of our second assumption and the inherent symmetries 
in the equations, all of the subsequent results can be reproduced if the starting formulation employs 
right-handed rotations instead. 

The first two components of r p / s can be found from the star-planet angular separation and a 
measurement of d, giving us the apparent separation (the distance of the planet from the star in 
the plane of the sky): 

" 1 
1 




(3) 



The third component of r p / s may also be observable if the detecting instrument can accurately 
measure the relative flux between the planet and its parent star, and some additional knowledge 
about the planet is assumed. Following iBrownl ((2005), we write the ratio of fluxes between planet 
and star as: 



Fr 
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with Amag = — 2.51og 10 Fr . 



(4) 



where p and R are the planet's albedo and radius, respectively, $ is the planet's phase function 
(|Sobolevlll975l ). and f3 is the star-planet-observer (phase) angle, given by: 
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However, since d> r, we can make the approximation: 
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(6) 



where z is the z axis component of r p / s . For the closest star to Earth (approximately 1.3 pc away), 
this approximation produces errors of less than 0.001%. As the direct detection completeness is a 
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function of only s and Fr, equations ([3]) and (HJ) are all that is needed to define the completeness. 
Our eventual goal is to find either the bivariate distribution function describing the direct detection 
completeness, or to find the marginal distribution functions of s and Fr, which can be sampled to 
generate the completeness distribution. To do so, we must first find expressions for the distributions 
of all of the parameters used in calculating these two quantities. 



3. Probability Density Functions of Orbital Parameters 



As described in £J2J the distributions of the Euler angles in the orbital parameter set are fixed by 
our assumption of uniform orbital orientation in space. The remaining two parameters, a and e, will 
be treated as unknowns, with probability density functions (PDFs) /g(e) and f a (a), respectively. 
These distribution s should ideally be generated via fits to observed data, which is an active area of 



current research. (jCumming et al.ll2008l : iBorucki et al. 



2010; 



Hogg et alJl20ld ) 



To find the distribution of planetary positions upon their orbits, we begin by noting that we 
can relate v to the eccentric anomaly (E) via the relationship: 
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which, in turn, is related to the mean anomaly (M) by the Kepler equation: 

M = E — esin(£) 



(7) 



(8) 



The mean anomaly is proportional to the area covered by r p / s , measured from periapsis passage, and 
is thus linear in time. A randomly timed observation therefore has equal probability of occurring 
at any value of M, so we let the mean anomaly be uniformly distributed in [0, 2n]. 

Let M, e and v be random variables with the distributions of mean anomaly, eccentricity and 
true anomaly, respe c tively , and define two functions: v = g(M,e) and M = h(i;,e). Following 



Larson and Shubertl (|1979l ). we can write the CDF for v as a marginalization of a conditional 



probability, 
F P (u) = P(g(M,e) < u) 



P(g(M,e) < u\e = e)Ue)de 



P(M<h(u,e))f B (e)de, (9) 



where the last step assumes independence between M and e, and that the orbits are closed so that 
e G [0, 1). We note that: 



P(M < h(u, e)) 
where /^(M) is the PDF of M: 

ffAM) = 



h(u,e) 



fu{M)&M 



else 



(10) 
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Therefore the PDF of v is: 



d 1 f 1 1% 1 f 1 (l — e 2 ) 

U(u) = ±F„(u) = — / °—h{e) de = - ^/ g (e) de , (12) 

dv 2vr 7 <9i/ 2?r 7 (1 + ecosz/ 



where /i is given by equations ((7j) and ([8]). 

Following the same procedure, but without substituting equation (|7|) into equation (|HJ), we find 
that the PDF of the eccentric anomaly is: 

fs(E) =J-^( E -esin(E)) f a (M)f s (e)de = i- J (1 - ecos£) / e -(e) de . (13) 
Similarly, if we rewrite the Kepler equation as: 

" = » -(^W 1 -!^ <14) 

and assume independence between a and e, we can write the PDF for the orbital radius as: 



/,!,•) 4ff , =/ g (e)de/ s (o)da. (15) 



7T 



Equations (fT2j) . (fT3j) . and (115]) are the most complete descriptions possible for the distributions of 
the orbital position and anomaly, without assuming specific distributions for the eccentricity and 
semi-major axis. In ^}5]we present examples where functions are selected for /g(e) and fa (a) and 
algebraic forms are calculated for fv{y) and ff(r). 



4. Probability Density Functions of Observed Quantities 

We derived the distribution functions of the true anomaly and orbital radius by assuming 
independence between the Keplerian orbital elements from which the anomaly and radius are 
calculated. However, the quantities observed by direct searches (i.e., s and Fr) are nonlinear 
functions of the Keplerian elements and r, so we cannot make the assumption of independence 
between their functional arguments for either of these. Fortunately, inspection of the phase angle 
reveals something interesting. Using the approximation in equation ([6]), we can write: 

(3 « cos^ 1 (cos 6 cos v — cos if) sin 6 sin v) (16) 

from which we define a new variable, 

x = cos /3 = cos 9 cos v — cos if> sin 6 sin v . (17) 

If it can be shown that x is uniformly distributed for all distributions of eccentricity, then this 
would prove that f3 is sinusoidally distributed regardless of the distribution of any other orbital 



parameter, given only our assumption of uniform orbital orientation. We do so by considering the 
characteristic function of x, 
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<p s (t)=M[e^)= / / ^Mu)f^)f B (9)d^d9du, (18) 



where E is the expectation value, and showing that it is equivalent to the characteristic function of 
a uniformly distributed variable. 

The assumption of uniform orbital orientation yields: 

m - { f i: ^ m 

m - ( r i: M 



so that the characteristic function of x becomes: 

1 

-oo JO JO 

Performing the integral over if;, 



oo r"K r2n 

<p s (t) = j- I / / e it( - cos0cosu - cos ^ sin0sin ^ dip sm6 d6f p (iy)dv (21) 



1 f'CO f'TT 

ips(t) = - / e i * CCBfloc " , J (tsmOwnv) sin0d0 f p (v) dv (22) 

2 J-oo Jo 

where Jq is the zeroth-order Bessel function of the first kin d. We n e xt pe rform the integral over 9 



using Gegenbauer's finite integral (see equation 12.14(1) in IWatsonl (119441 )): 



J\ itcoBeamv J b _i (tsm0smv)C b a (cos9)sm b+ ^ 0d9 = y^ a s i n b -§ „ c*(cos v) J a+b {t) (23) 

where C b is a Gegenbauer polynomial. Choosing a = and b = | yields C b (x) = 1 and the integral 
becomes: 

J\ ltcos9cosu J (tsm9smv)sm6d9 = ] J^Ji(t) (24) 



'o 

Since the half-order Bessel function is defined as: 



/ 2 

Ji(t) = \—smt (25) 

2 \ lit 



the characteristic function becomes: 
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By definition, 



and thus, 



j — ( 
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(28) 
(29) 



If we define a random variable y ~ £/([—!, 1]), its characteristic function will be: 
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e lty fy{y)dy 
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it 
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(30) 



As equations (|29p and (|30p are identical, x and y have the same characteristic function, and thus 
must have the same underlying distribution. Since x is uniformly distributed on [—1,1], j3 is 
sinusoidally distributed in [0, tt) for any population of uniformly oriented orbits, independent of the 
specific distribution of e and a. This finding greatly simplifies our derivation of the distribution 
functions of the flux ratio and apparent separation. 

Returning now to the definition of flux ratio in equation (JH), and assuming that we are in- 
terested in a particular planet type such that pR 2 is approximately constant, we see that the flux 
ratio is a function of the product of two variables, m = 3>(/3) and n = r~ 2 . If we define a random 
variable k = mn, then it can be shown that: 

coo ^ 
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(31) 



if m and n are independent. ( Larson and Shubertlll979 ) As we have demonstrated the independ- 
of P from f, equation (i3~TT) is applicable. Using equation (fl5l) . the CDF of n is given by: 
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The PDF of n is thus 
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where C = a 2 (e 2 - l)n. The PDF of m is dependent on the exact form of but can also be 
evaluated analytically when $ is invertible. In those cases, 



frn(m) = f-^- L (m)) 



-^-<E> 1 (m) 
am 



±sin ($ _1 (m)) |a^$ -1 ("i)| $ _1 (m) G [0,7r) 







else 



(34) 



This formulation is still problematic when dealing with transcen dental functions such as the 
commonly used phase function of a Lambertian sphere (|Sobolevlll975l ). 



tt$(/3) = sin /3 + (vr - 0) cos j3 . 



(35) 



but, since the domain of the function is restricted to f3 G [0, tt], the range of $ is single- valued 
(G [1,0]), and the function is invertible. Expanding about vr/2, equation (I35p becomes: 



, 1 ^ / vrxfc 4 /2(A;-1)\ 2d fe 

^ = - + S a H /3 " 9) where «fc 



fc=i 



2kl 



TT 



(36) 



and di : 3 = —1,1,1, with dh = dk--\ dk-?d k-a for > 3, such that you get the series d{ = 
{— 1, 1, 1, — 1, — 1, 1, 1, — 1, — 1 . . .}. Following iMorse and Feshbachl (|1953l ). we can express the in- 
verse function as the series: 
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and the space X of sets x is defined as: 
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(38) 



fc-l 



1= i£ N fc_1 : ^ ixi = A; — 1 > with |x|=5^: 



(39) 



The inverse series converges to machine double-precision in a finite number of terms, except at the 
endpoints (0 and 1), which themselves do not need to be computed as they map exactly to f3 = tt 
and /3 = 0, respectively. 

Since Fr = pR 2 k, equations ([3T]) - ([39]) allow us to write: 



1 



pR 2 J_ 



fn(n) 



■ cos 



II 



dn , (40) 



which represents the PDF of the planet-star flux ratio for a population of planets with constant 
pR 2 values (i.e., Earth- twins) modeled as Lambertian reflectors. 
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Finally, we consider the apparent separation. The same approximation used in equation ([6]) 
allows us to write: 

s = r sin j3 . (41) 



If we let I = sin/3, following equation (|34[) : 
/K/) = /^(sin- 1 (/)) 



Asin-(/) 



I G [0,1) 
else 



(42) 



Returning to equations f)15|) and (|3ip and now letting m = r and n = I, we have: 
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(43) 



which represents the probability density function of the apparent separation. We now have ex- 
pressions for the distribution functions of both of the observed quantities of direct imaging planet 
searches. While the expressions are quite complex, they are numerically integrable, and can be 
used to check and improve the efficiency of the sampling of the completeness function. With an 
assumed /g(e) and fa (a), equations PU|) and (|4"3j) . via marginalization and joint sampling, can 
directly produce the completeness distribution. 



5. Validation of Derived Distributions 

As a check on the derived distribution functions for the true anomaly and orbital radius, we 
first consider the simplest possible case, where both eccentricity and semi-major axis are uniformly 
distributed: 

£ I \ j (o>max ^mira) Q G l&max i Q"min] / . r \ 

Ua) = else (45) 



Equation (112]) then becomes: 



1 rl (l-e 2 ) 1 
2vr J (1 + ecosz/) 



- 2 Fi(l,2;7/2,cos 2 i/)— — j- (4 sin v + cos 2u - 3) (46) 

07T lb cos 4 V 
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where 2-P1 is the Gaussian hypergeometric function. Equation (|15p becomes: 
1 



Mr) 



7ra (o 

r, 



[2r log (Ci + y/r) + a log (a — r 
+ ia log 



+ ia log 



/ 4VF(2C 2 + C 1 +iV?) 
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Vr(-2C7 2 + C7 1 + iVr) 



2r log (C2 — \/r — a) 
-2alog(-2(r + VFCi)) (47) 
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where C\ = \/2a — r and C 2 = 



r. 



To check these equations, we generat e one million IIP samples each from the uniform distri- 
butions in equations (|4"4"|) . (j4"5j) and ([TT]) , (jPress et al.lll992l ) From these, we calculate the eccentric 
anomalies via Newton- Raphson iteration applied to equation ([8]), and then calculate v and r via 
equations ([7]) and ([2]) , respectively. The sample PDFs are calculated for these two parameters and 
compared with the results of our closed- form solutions, as shown in Figures 1(a) and 1(b). These 
plots demonstrate excellent agreement, validating the equations. 

While uniform distributions are a useful check, th e semi-major axes of currently found planets 
appear more likely to be logarithmically distributed. (jCurrid 120091 ) As an additional test, we will 
assume that a is uniformly distributed in log 10 (a) for a G [10 Cmi ™, 10 Cma:r ] so that its PDF is: 



Ha) 



(Aclog(lO)a) 




-1 



a G [10 c "»™,10 c 
else 



(48) 



where Ac = c max — c m i n . If we retain the same uniform distribution for eccentricity, the true 
anomaly distribution remains the same, but the orbital radius density becomes: 

'32r 3 , 
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2vrlog(10)Aca 2 r 



aC\\[r — a log (a 



2ia log 



CM 



+2a 2 log -4 C\r 2 + r 2 + 2r 2 log 



C 2 ~ Vr~^ 



io c 



(49) 



where C\ and C 2 are defined as in equation (|47p . We repeat the simulation, this time generating our 
sample of semi-major axes by exponentiating 10 with a sample of one million uniformly distributed 
random values between -1 and 1 (thereby setting Ac to 2), and re-calculating the PDF for / f . This 
is compared with the results from equation (|49p in Figure 1(c), and shows excellent agreement. 

We can also empirically check the assertion that the distribution of (3 is independent of the 
distribution of e. We now generate two sets of IID samples, one using the uniform distribution of 
e, and one distributed via a step function, 



Me) 



1.3 e G [0,0.7) 
0.3 e G [0.7,1] 
else 



(50) 
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where accurate sampling is achieved via simple rejection sampling. The two sample sets are used to 
generate two different sample distributions of (3. Figure [2] shows Q-Q plots comparing the quantiles 
of these two samples to the theoretical quantiles of a sinusoidal distribution. Quantiles for the ideal 
sinusoidal distribution are calculated by evaluating its inverted CDF at regularly spaced intervals, 
while quantiles of the two simulated distributions are calculated by ordering and binning the gen- 
erated values. As the resulting Q-Q plots follow the 45° line, this demo nstrates that the simulated 



distr ibutions are identical to the theoretical sinusoidal distribution. (IGibbons and Chakraborti 
2003h 



6. Conclusions 



Starting with a uniform distribution of orbital orientations and arbitrary PDFs for eccentricity 
and semi-major axis, we have derived expressions for the distributions of the remaining Keplerian 
orbital elements, and of parameters used to describe direct exoplanet detections. We have demon- 
strated the independence of phase angle from the distributions of the Keplerian orbital elements, 
allowing for the calculation of the distributions of both apparent separation and flux ratio, which 
can be directly sampled to generate the completeness distribution. At the same time, we have 
provided fully analytic forms for the distributions of Keplerian orbital elements based on uniform 
distributions of eccentricity and semi-major axis, as well as logarithmically distributed semi-major 
axes, assumptions often made in research related to planet-finding mission planning and analysis. 
These forms should allow researchers to increase the efficiency of their simulations, and to empir- 
ically check for global errors generated by under-sampling. The same equations and procedures 
described here will also be useful for other statistical work associated with planet-finding, including 
calculating the likelihood of transits or inferring the true distribution of orbital parameters from 
doppler spectroscopy surveys. 



The authors gratefully acknowledge the input of Robert Vanderbei, whose comments and 
suggestions were invaluable, and Robert Brown, whose work has so much informed our own. 



REFERENCES 

W. J. Borucki, D. Koch, G. Basri, N. M. Batalha, T. Brown, D. A. Caldwell, J. Caldwell, 
J. Christensen-Dalsgaard, W. Cochran, E. DeVore, E. Dunham, A. Dupree, T. Gautier, 
J. Geary, R. Gilliland, A. Gould, S. Howell, J. Jenkins, H. Kjeldsen, Y. Kondo, D. Latham, 
J. Lissauer, G. Marcy, S. Meibom, D. Monet, D. Morrison, D. Sasselov, and J. Tarter. Ke- 
pler Planet Detection Mission: Introduction and First Results. In American Astronomical 
Society Meeting Abstracts, volume 215 of American Astronomical Society Meeting Abstracts, 
page 101.01, January 2010. 



-13- 



R. A. Brown. Obscurational completeness. ApJ, 607:1003-1017, 2004. 

R. A. Brown. Single-visit photometric and obscurational completeness. ApJ, 624:1010-1024, 2005. 

R. A. Brown. On the completeness of reflex astrometry on extrasolar planets near the sensitivity 
limit. The Astrophysical Journal, 699(1):711-715, 2009. 

R.A. Brown and R. Soummer. New completeness methods for estimating exoplanet discoveries by 
direct detection. The Astrophysical Journal, 715:122, 2010. 

A. Cumming, R. P. Butler, G. W. Marcy, S. S. Vogt, J. T. Wright, and D. A Fischer. The 
Keck planet search: detectability and the minimum mass and orbital period distribution of 
extrasolar planets. PASP, 120:531-554, 2008. 

T. Currie. On the semimajor axis distribution of extrasolar gas giant planets: Why hot jupiters 
are rare around high- mass stars. The Astrophysical Journal Letters, 694:L171, 2009. 

J.D. Gibbons and S. Chakraborti. Nonparametric statistical inference. CRC, 2003. 

A. Gould, S. Dong, BS Gaudi, A. Udalski, IA Bond, J. Greenhill, RA Street, M. Dominik, T. Sumi, 
MK Szymahski, et al. Frequency of solar-like systems and of ice and gas giants beyond 
the snow line from high-magnification microlensing events in 2005-2008. The Astrophysical 
Journal, 720:1073, 2010. 

W.K. Hastings. Monte Carlo sampling methods using Markov chains and their applications. 
Biometrika, 57(1):97-109, 1970. 

D.W. Hogg, A.D. Myers, and J. Bovy. Inferring the eccentricity distribution. Arxiv preprint 
arXiv:1008.4146, 2010. 

H. J. Larson and B. O. Shubert. Probabilistic Models in Engineering Sciences Vol. 1. John Wiley 
& Sons, New York, 1979. 

P. M. Morse and H. Feshbach. Methods of Theoretical Physics. McGraw-Hill New York, 1953. 

W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical recipes in C: The 
Art of Scientific Computing. Cambridge University Press, 1992. 

D. Savransky, N. J. Kasdin, and E. Cady. Analyzing the designs of planet finding missions. PASP, 
122:401-419, April 2010. 

V. V. Sobolev. Light Scattering in Planetary Atmospheres. Pergamon Press, New York, 1975. 

J. P. Vinti. Orbital and Celestial Mechanics. Charles Stark Draper Laboratory, Inc., Cambridge, 
Mass., 1998. 

G. N. Watson. A treatise on the theory of Bessel functions. Cambridge University Press, 1944. 



-14- 




A 



1 2 3 4 5 6 0.5 1 1.5 2 2.5 3 




r 



(c) /r, log a 

Fig. 1. — Comparison of empirical (black line) and derived (gray line) distributions assuming a 
uniform distribution for e and uniform and log distributions for a. For uniformly distributed semi- 
major axis we assumed a £ [0.5, 1.5], leading to r G (0, 3], and for a uniform in log a, we assumed 
a G [0.1, 10] leading tore (0, 20]. 
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